/*
 * calculate_gradient.cpp
 *
 *  Created on: 13/09/2012
 *      Author: rogsoares
 */


#include "EBFV1_elliptic.h"

namespace PRS{
	double EBFV1_elliptic::pressureGradient(TMesh* pTMesh_data){
		int dim = pTMesh_data->getMeshDim();
		int k, nrows, ncols, dom;
		for (dom=0; dom<pMatVec->ndom; dom++){

			// reset gradient data grad_p = 0 for domain 'dom'
			pTMesh_data->vertices_list[dom].getSize(nrows,ncols);
			for (int j=0; j<nrows; j++){
				VertexData* vdata = pTMesh_data->vertices_list[dom].getValue(j,0);
				for(int i=0;i<dim;i++){
					vdata->grad_p[i] = .0;
				}
			}

			// calculate gradient: all edges
			pTMesh_data->edge_list[dom].getSize(nrows,ncols);
			for (int i=0; i<nrows; i++){
				EdgeData* edata = pTMesh_data->edge_list[dom].getValue(i,0);
				VertexData* vertex_0 = edata->vertex_0;
				VertexData* vertex_1 = edata->vertex_1;
				double val = 0.5*(vertex_0->p + vertex_1->p);
				for (k=0; k<dim; k++){
					vertex_0->grad_p[k] += val*edata->Cij[k];
					vertex_1->grad_p[k] += -val*edata->Cij[k];
				}
			}

			// calculate gradient: only boundary elements (2-D: edges or 3-D: triangles
			pTMesh_data->bdry_elem_list[dom].getSize(nrows,ncols);
			for (int i=0; i<nrows; i++){
				BdryElementData* edata = pTMesh_data->bdry_elem_list[dom].getValue(i,0);
				VertexData* vertex_0 = edata->vertex_0;
				VertexData* vertex_1 = edata->vertex_1;
				double val = 0.5*(vertex_0->p + vertex_1->p);
				for (k=0; k<dim; k++){
					vertex_0->grad_p[k] += ((5.*vertex_0->p + vertex_1->p)/6.0)*edata->Dij[k];
					vertex_1->grad_p[k] += ((vertex_0->p + 5.*vertex_1->p)/6.0)*edata->Dij[k];
				}
			}

			// calculate gradient: weighting by node volume
			pTMesh_data->vertices_list[dom].getSize(nrows,ncols);
			for (int i=0; i<nrows; i++){
				VertexData* vdata = pTMesh_data->vertices_list[dom].getValue(i,0);
				for (int k=0; k<dim; k++) {
					vdata->grad_p[k] /= vdata->volume;
				}
			}
		}
		return 0;
	}
}
